/*-------------------------------------------------------------------------*/
/**
  @file		grid_func.c
  @author	Renaud Wanchoor
  @update date	March 2007 by Xiaochen Wu
  @brief	grid utilities

  Management of the structure that store the subdomain connectivity graph,
  provides also the swaaping stage and functions for basic grid management:
  - Initialization
  - Finalization
  - read from a file

*/
/*--------------------------------------------------------------------------*/

#include <plib.h>
#include "pagrid.h"

#define LW_CEIL 1000

/*********************************************/
/*        swapping processors  tools         */
/*********************************************/
#define XADJ(npro,x,y) (x*((npro<<1)-x-3)-2+(y<<1))>>1 /**< \brief compute graph->adjidx for a fully connected graph, work only if x<y */
#define XADJT(npro,x,y) x<y?XADJ(npro,x,y):XADJ(npro,y,x) /**< \brief compute graph->adjidx for a fully connected graph */


static double swapGain  (SwapCtrl*,int prov,int prow);
static void   swapUpdateGraph(SwapCtrl* ctrl,GraphType* graph,int edge);
static int    swapUpdate(SwapCtrl*,int edge);

/*********************************************/
/*            connectivity tools             */
/*********************************************/
static void UpdateCon(GraphType* graph,int from,int to,int ewgt);


/****************************************************************************/
/****************************************************************************/
/**                                                                        **/
/**                      FUNCTIONS (connectivity)                          **/
/**                                                                        **/
/****************************************************************************/
/****************************************************************************/
/* graph->app contains the graph that represents the subdomain connectivity */

/**
 * \brief initialize the subdomain connectivity graph
 *
 * copy graph to GridType::app the graph of the subdomain connectivity, graph needs
 * to be the coarsest graph
 *
 * \params grid computational grid structure
 * \params graph coarsest graph structure
 */
void IniProcessorCon(GridType* grid,GraphType* graph){
  
  int nproc      = grid->npro;
  int nnedges    = SQUARE(nproc);
  GraphType* con = grid->app;
  Mem* mem       = &(grid->app->mem);

  int *nadjidx,*nadjwgt,*nadjncy,*nvwgt,*npro;

  int  gnedges   = graph->nedges;
  int* gadjidx   = graph->adjidx;
  int* gadjncy   = graph->adjncy;
  int* gadjwgt   = graph->adjwgt;
  int* gvwgt     = graph->vwgt;
  int* gpro      = graph->pro;
  
  con->nvtxs  = graph->nvtxs;
  con->nedges = graph->nedges;

  MEM_Need(mem,3*nproc+2*nnedges+1,int);
  MEM_Allocate(mem);
      
  con->vwgt   = nvwgt   = MEM_Request(mem,nproc    , int);
  con->adjidx = nadjidx = MEM_Request(mem,nproc + 1, int);
  con->adjncy = nadjncy = MEM_Request(mem,nnedges  , int);
  con->adjwgt = nadjwgt = MEM_Request(mem,nnedges  , int);
  con->pro    = npro    = MEM_Request(mem,nproc    , int);

  memcpy(nvwgt,  gvwgt,  nproc    *sizeof(int));
  memcpy(nadjidx,gadjidx,(nproc+1)*sizeof(int));
  memcpy(nadjncy,gadjncy,gnedges  *sizeof(int));
  memcpy(nadjwgt,gadjwgt,gnedges  *sizeof(int));
  memcpy(npro,   gpro,   nproc    *sizeof(int));
}

/**
 * \brief update an edge in the subdomain connectivity graph
 *
 * update the edge weight between from and to, if the edge was not existing it
 * is created, if after the update the edge weight between from and to is 0 it
 * is removed
 *
 * \param graph graph to be updated
 * \param from endpoint 1
 * \param to endpoint 2
 * \param ewgt edge weight to add between from and to
 */
void UpdateCon(GraphType* graph,int from,int to,int ewgt){
  int i,e = -1,istart,iend;

  int  nvtxs  = graph->nvtxs; 
  int  nedges = graph->nedges;
  int* adjidx = graph->adjidx;
  int* adjncy = graph->adjncy;
  int* adjwgt = graph->adjwgt;

  /* search existing link */
  istart = adjidx[from];
  iend   = adjidx[from+1];

  for(i=istart;i<iend;i++){
    if(adjncy[i] == to){
      e = i; 
      break;
    }
  }

  /* create link */
  if(e == -1){
    for(i=from;i<nvtxs;i++)
      adjidx[i+1] += 1;
    for(i=nedges;i>iend;i--){
      adjncy[i] = adjncy[i-1];
      adjwgt[i] = adjwgt[i-1];
    }
    adjncy[iend] = to;
    adjwgt[iend] = ewgt;
    graph->nedges++;

    return;
  }

  /* update link and remove it if is 0 */
  if((adjwgt[e] += ewgt) == 0){
    for(i=from;i<nvtxs;i++)
      adjidx[i+1] -= 1;
    for(i=e;i<nedges-1;i++){
      adjncy[i] = adjncy[i+1];
      adjwgt[i] = adjwgt[i+1];
    }
    graph->nedges--;
  } 
}

/**
 * \brief update the subdomain connectivity graph when moving vertex vtx
 *
 * Update the subdomain connectivity (grid->app) when vertex vtx is moved from
 * processor from to processor to
 *
 * \param grid computational grid
 * \param graph application graph
 * \param vtx vertex moved
 * \param from origin of the vertex
 * \param to destination of the vertex
 */
void UpdateProcessorCon(GridType* grid,GraphType* graph,
			int vtx,int from, int to){
  int i,istart,iend;
  int w,prow,vprow;
  
  int* gadjidx = graph->adjidx;
  int* gadjncy = graph->adjncy;
  int* gadjwgt = graph->adjwgt;
  int* gvwgt   = graph->vwgt;
  int* gpro    = graph->pro;

  int* nvwgt   = grid->app->vwgt;
  int* rpro    = grid->rpro;

  int vfrom    = rpro[from];
  int vto      = rpro[to];

  nvwgt[vfrom] -= gvwgt[vtx];
  nvwgt[vto]   += gvwgt[vtx];  
  
  istart = gadjidx[vtx];
  iend   = gadjidx[vtx+1];
  for(i=istart;i<iend;i++){
    w    = gadjncy[i];
    prow = gpro[w];
    
    if(prow == from){
      // add to processor connectivity
      UpdateCon(grid->app,vto,vfrom,gadjwgt[i]);
      UpdateCon(grid->app,vfrom,vto,gadjwgt[i]);
      continue;
    }
    
    if(prow == to){
      // remove from processor connectivity
      UpdateCon(grid->app,vto,vfrom,-gadjwgt[i]);
      UpdateCon(grid->app,vfrom,vto,-gadjwgt[i]);
      continue;
    }
    
    // update connectivity
    vprow = rpro[prow];
    UpdateCon(grid->app,vfrom,vprow,-gadjwgt[i]);
    UpdateCon(grid->app,vprow,vfrom,-gadjwgt[i]);
    UpdateCon(grid->app,vto,vprow,gadjwgt[i]);
    UpdateCon(grid->app,vprow,vto,gadjwgt[i]);
  }   
}

/****************************************************************************/
/****************************************************************************/
/**                                                                        **/
/**                      FUNCTIONS (swapping processors)                   **/
/**                                                                        **/
/****************************************************************************/
/****************************************************************************/
/**
 * \brief swapping stage
 *
 * swap processors assignment of current partition that decrease the total
 * estimated execution time without increasing the maximum execution time
 *
 * \param grid coputational grid
 * \param graph application graph
 */
void SwapProcessors(GridType* grid,GraphType* graph){
  int i,j,e;
  int prov,prow,edge;

  Mem       mem;
  SwapCtrl  ictrl;
  SwapEdge* edges;

  int       npro = grid->npro;
  double*   exec = grid->exec;

  int     nedges = npro*(npro-1)/2;
  
  double gain;

  Heap h_gain;			/* gain heap             */

  /**********************************/
  /* Init                           */
  /**********************************/
  /* prepare */
  MEM_Init(&mem);
  HeapInit(&mem,&h_gain,nedges);
  MEM_Need(&mem, npro,   double);
  MEM_Need(&mem, nedges, SwapEdge);
  
  /* allocate */
  MEM_Allocate(&mem);

  ictrl.pexec         = MEM_Request(&mem, npro,   double);
  ictrl.edges = edges = MEM_Request(&mem, nedges, SwapEdge);

  /* initialize */
  ictrl.heap = &h_gain;
  ictrl.g    = grid->app;
  ictrl.n    = grid;

  e=0;
  for(i=0;i<npro;i++)for(j=i+1;j<npro;j++){
    edges[e].prov = i;
    edges[e].prow = j;
    edges[e].pos  = NULL;
    e++;
  }

  ictrl.tmax = 0;
  for(i=0;i<npro;i++){
    if(exec[i]>ictrl.tmax)
      ictrl.tmax = exec[i];
  }

  /************************************/
  /* Initial state                    */
  /************************************/
  /* compute gain when switching processors assignement */
  for(prov=0;prov<npro;prov++){
    for(prow=prov+1;prow<npro;prow++){
      gain = swapGain(&ictrl,prov,prow);
      
      HeapInsert(&h_gain, gain);
    }
  }
  
  /****************************************/
  /* Refinement                           */
  /****************************************/
  while(HeapSeeHead(&h_gain) > 0){

    edge = HeapRemoveHead(&h_gain);
    if(swapUpdate(&ictrl,edge)){
      swapUpdateGrid(&ictrl, graph, grid, edge);
      swapUpdateGraph(&ictrl,graph,edge);
    }
  }
  
  MEM_Free(&mem);
}

/**
 * \brief compute the gain of when swapping two processors
 *
 * compute the gain of switching prov with prow and put
 * update in execution time of each processor into ctrl->pexec
 *
 * \param ctrl swapping ctrl structure
 * \param prov processor to be swapped
 * \param prow processor to be swapped
 */
double swapGain(SwapCtrl* ctrl,int prov,int prow){
  int x,j;
  int prox;

  int *adjidx      = ctrl->g->adjidx;
  int *adjncy      = ctrl->g->adjncy;
  int *adjwgt      = ctrl->g->adjwgt;
  int *vwgt        = ctrl->g->vwgt;
  int *pro         = ctrl->g->pro;
  
  double **ratio   = ctrl->n->ratio;
  double *exec	   = ctrl->n->exec;	/* overall execution time */
  int npro	   = ctrl->n->npro;
  double* prospeed = ctrl->n->prospeed;

  int **latency    = ctrl->n->latency;

  int *rpro        = ctrl->n->rpro;
  double *pexec    = ctrl->pexec;	/* update to do on execution time */

  int v		   = rpro[prov];
  int w		   = rpro[prow];
  double gain      = 0;

  /******************************/
  /* gain info already existing */
  /******************************/
  for(j=0;j<npro;j++)
    pexec[j] = 0;

  /******************************************************/
  /* gain assigning prow to v and adjacent consequences */
  /******************************************************/
  pexec[prow] = exec[prow] - vwgt[v]*prospeed[prow];
  for(j=adjidx[v];j<adjidx[v+1];j++){
    x = adjncy[j];

    if(x == w)		/* processors were swapped */
      prox = prov;
    else
      prox = pro[x];
    
    pexec[prow] -= adjwgt[j]*ratio[prow][prox];
    pexec[prow] -= latency[prow][prox];

    if(x != w){
      pexec[prox] += adjwgt[j]* (ratio[prox][prov] - ratio[prox][prow]);
      pexec[prox] += latency[prox][prov] - latency[prox][prow];
    }
  }

  /******************************************************/
  /* gain assigning prow to v and adjacent consequences */
  /******************************************************/
  pexec[prov]    = exec[prov] - vwgt[w]*prospeed[prov];

  for(j=adjidx[w];j<adjidx[w+1];j++){
    x = adjncy[j];

    if(x == v)		/* processors were swapped */
      prox = prow;
    else
      prox = pro[x];
    
    pexec[prov]   -= adjwgt[j]*ratio[prov][prox];
    pexec[prov]   -= latency[prov][prox];

    if(x != v){
      pexec[prox] += adjwgt[j]* (ratio[prox][prow] - ratio[prox][prov]);
      pexec[prox] += latency[prox][prow] - latency[prox][prov];
    }
  }

  /**************/
  /* total gain */
  /**************/
  for(j=0;j<npro;j++)
    gain += pexec[j];

  return gain;
}
/**
 * \brief update the graph when two processors are swapped
 *
 * update the application graph when the two processors that are
 * endpoints of the edge are swapped
 *
 * \param ctrl swapping ctrl structure
 * \param graph application graph
 * \param edge edge id in a fully connected graph
 */
void swapUpdateGraph(SwapCtrl* ctrl,GraphType* graph,int edge){
  int i;

  int nvtxs       = graph->nvtxs;
  int prov	  = ctrl->edges[edge].prov;
  int prow	  = ctrl->edges[edge].prow;
  int* pro        = graph->pro;

  for(i=0;i<nvtxs;i++){
    if(pro[i] == prov){
      pro[i] = prow;
      continue;
    }
    if(pro[i] == prow){
      pro[i] = prov;
    }
  }
}

void swapUpdateGrid(SwapCtrl* ctrl, GraphType* graph, GridType* grid, int edge){
  int v, i, istart, iend;
  int x, prox;

  int* adjidx     = graph->adjidx;
  int* adjncy     = graph->adjncy;
  int* pro        = graph->pro;
  int nvtxs       = graph->nvtxs;
  int npro        = grid->npro;
  int** connect   = grid->connect;
  int **convtxs   = grid->convtxs;
  int  *conproc   = grid->conproc;
  int prov        = ctrl->edges[edge].prov;
  int prow        = ctrl->edges[edge].prow;

  /* Updage connecting information of the grid */
  for(v = 0; v < nvtxs; v++){
    if(pro[v] == prov){
      istart = adjidx[v];
      iend = adjidx[v+1];
      for(i = istart; i < iend; i++){
        x = adjncy[i];
        prox = pro[x];
        if(prox != prov && prox != prow){
          connect[prov][prox]--;
          connect[prox][prov]--;
          connect[prow][prox]++;
          connect[prox][prow]++;
        }
      }
    }
    if(pro[v] == prow){
      istart = adjidx[v];
      iend = adjidx[v+1];
      for(i = istart; i < iend; i++){
        x = adjncy[i];
        prox = pro[x];
        if(prox != prow && prox != prov){
          connect[prow][prox]--;
          connect[prox][prow]--;
          connect[prov][prox]++;
          connect[prox][prov]++;
        }
      }
    }
  }
  for(i = 0; i < npro; i++){
    conproc[i] = convtxs[prov][i];
    convtxs[prov][i] = convtxs[prow][i];
    convtxs[prow][i] = conproc[i];
  }
  for(i = 0; i < npro; i++){
    conproc[i] = convtxs[i][prov];
    convtxs[i][prov] = convtxs[i][prow];
    convtxs[i][prow] = conproc[i];
  }
}

/**
 * \brief update the swap heap when two processors are swapped
 *
 * update the heap when the two processors that are
 * endpoints of the edge are swapped
 *
 * \param ctrl swapping ctrl structure
 * \param edge edge id in a fully connected graph
 */
int swapUpdate(SwapCtrl* ctrl,int edge){
  int i,j,k,x,z,prox,proz;
//  int istart, iend;

  double gain;
  double old_tmax;
  Heap* h	  = ctrl->heap;

  int* adjidx     = ctrl->g->adjidx;
  int* adjncy     = ctrl->g->adjncy;
  int* adjwgt     = ctrl->g->adjwgt;
  int* pro	  = ctrl->g->pro;

  int* rpro	  = ctrl->n->rpro;
  int npro	  = ctrl->n->npro;
  double** ratio  = ctrl->n->ratio;
  double* exec    = ctrl->n->exec;
  int **latency   = ctrl->n->latency;
//  int **connect   = ctrl->n->connect;

  int prov	  = ctrl->edges[edge].prov;
  int prow	  = ctrl->edges[edge].prow;
  int v 	  = rpro[prov];
  int w 	  = rpro[prow];

  /* update exec time and invert gain of switched edge */
  /* filter swap that increase tmax */
  swapGain(ctrl,prov,prow); /* this update ctrl->pexec */
  for(i=0;i<npro;i++){
    exec[i] -= ctrl->pexec[i];
  }

  old_tmax   = ctrl->tmax;
  ctrl->tmax = 0;
  for(i=0;i<npro;i++){
    if(exec[i] >= ctrl->tmax)
      ctrl->tmax = exec[i];
  }

  if(old_tmax < ctrl->tmax){
    //NOTICE((1,"swap increase tmax (%10.2f=>%10.2f)\n",old_tmax,ctrl->tmax));
    for(i=0;i<npro;i++){
      exec[i] += ctrl->pexec[i];
    }
    ctrl->tmax = old_tmax;
    return 0;
  }
 
  /* switch processors keeping old assignement in prov and prow */
  pro[v] = prow;
  pro[w] = prov;
  rpro[prov] = w;
  rpro[prow] = v;

  /**************************************************/
  /* Update gains in the heap "it's all in the hip" */
  /**************************************************/
  /* consequences of switching prov with prow       */
  /**************************************************/
  /* update all (prox;prov) couples exept when prox = prow (already done) */
  /* this also update edges that has an end common with w                 */
  for(prox=0;prox<prov;prox++){
    if(prox == prow)
      continue;
    j = XADJ(npro,prox,prov);
    
    /* recompute gain (need to go throw all edges of x and w anyway) */
    gain = swapGain(ctrl,prox,prov);

    /* update heap */
    HeapUpdate(h, j, gain);
  }
  
  for(prox = prov+1;prox<npro;prox++){
    if(prox == prow)
      continue;
    j = XADJ(npro,prov,prox);

    gain = swapGain(ctrl,prox,prov);
    HeapUpdate(h, j, gain);
  }

  /* update edge gain for those who as a link (i) with v */
  for(i=adjidx[v];i<adjidx[v+1];i++){
    x = adjncy[i];
    if(x == w)			/* direct link with w already done */
      continue;
    prox = pro[x];

    for(proz=0;proz<npro;proz++){
      if(proz == prov || proz == prow || proz == prox)
	continue;
      z = rpro[proz];
      
      k = XADJT(npro,prox,proz);
      gain  = HeapGetGain(h,k);
      /* update gain other the link (remove prov add prow) */
      gain  += adjwgt[i]*(- ratio[prox][prov] + ratio[prox][prow] 
			  + ratio[proz][prov] - ratio[proz][prow]
			  - ratio[prov][prox] + ratio[prov][proz]
			  + ratio[prow][prox] - ratio[prow][proz])
                          - latency[prox][prov] + latency[prox][prow]
                          + latency[proz][prov] - latency[proz][prow]
                          - latency[prov][prox] + latency[prov][proz]
                          + latency[prow][prox] - latency[prow][proz];
      HeapUpdate(h, k, gain);
    }
  }

  /**************************************************/
  /* consequences of switching prow with prov       */
  /**************************************************/
  /* update all (prox;prow) couples exept when prox = prov (already done) */
  /* this also update edges that has an end common with v                 */
  for(prox=0;prox<prow;prox++){
    if(prox == prov)
      continue;
    
    j = XADJ(npro,prox,prow);
    
    /* recompute gain (need to go throw all edges of x and v anyway) */
    gain = swapGain(ctrl,prox,prow);

    /* update heap */
    HeapUpdate(h, j, gain);
  }
  
  for(prox = prow+1;prox<npro;prox++){
    if(prox == prov)
      continue;

    j = XADJ(npro,prow,prox);

    gain = swapGain(ctrl,prox,prow);
    HeapUpdate(h, j, gain);
  }
  
  /* update edge gain for those who as a link (i) with v */
  for(i=adjidx[w];i<adjidx[w+1];i++){
    x = adjncy[i];
    if(x == v)				/* direct link with w already done */
      continue;
    prox = pro[x];
    
    for(proz=0;proz<npro;proz++){
      if(proz == prov || proz == prow || proz == prox)
	continue;
      
      z    = rpro[proz];
      k    = XADJT(npro,proz,prox);
      gain = HeapGetGain(h,k);

      /* update gain other the link (remove prow add prov) */
      gain += adjwgt[i]*(- ratio[prox][prow] + ratio[prox][prov]
			 + ratio[proz][prow] - ratio[proz][prov]
			 - ratio[prow][prox] + ratio[prow][proz]
			 + ratio[prov][prox] - ratio[prov][proz])
                         - latency[prox][prov] + latency[prox][prow]
                         + latency[proz][prov] - latency[proz][prow]
                         - latency[prov][prox] + latency[prov][proz]
                         + latency[prow][prox] - latency[prow][proz];
      HeapUpdate(h, k, gain);
    }
  } 

  return 1;
}


/****************************************************************************/
/****************************************************************************/
/**                                                                        **/
/**                      FUNCTIONS (general grid functions)                **/
/**                                                                        **/
/****************************************************************************/
/****************************************************************************/
/**
 * \brief Initialize the grid
 *
 * dynamically allocate grid (need to be freed with LIBN_Free()) and
 * initialize the structure
 */
GridType* LIBN_Init(void){
  GridType *g;

  g = U_malloc(1,GridType,"LIBN_Init code 1\n");
  g->netcost = g->realnetcost = NULL;
  g->penalty = NULL;
  g->proratio = g->proportion = NULL;
  g->connect = NULL;
  g->maxdegree = g->npro = 0;  
  g->app = NULL;
  g->rpro = NULL;  

  MEM_Init(&(g->mem));

  return g;

}

/**
 * \brief read a grid file
 *
 * dynamically allocate grid (need to be freed with LIBN_Free) and
 * initialize the structure with the data contained in the file
 *
 * \param fp file pointer to the file that describes the grid
 */
GridType* LIBN_Read(FILE* fp, int NBytes){
  int i,j,k,x;
  int a, b;
  int dist;
  int ltcy;
  int istart,iend;

  char line[MAXLINE];
  char *tok1,*tok2,*tok3;
  int linenb = 1;
  
  GridType* g;
  
  int npro = 0;
  int nedges = 0;  			
  int **realnetcost,**netcost;
  int **connect;
  int **latency;
  int **latency1;
  int **convtxs;
  int **convtxs2;
  int *conproc;
  double **ratio,*proportion;
  double *prospeed;
  int maxproportion = 0;
  int *adjidx,*adjncy,*adjwgt,*adjlcy;
  Mem *mem;
  int max_latency = 0;
  int range = 1;

  /* reference ratio readen from file */
  int r_from,r_to;
  double r_val;
  double r_ref;  
  
  g = LIBN_Init();
  mem = &(g->mem);
  
  /******************************************/
  /* read two first lines                   */
  /******************************************/
  /* read first line (number of processors and edges number) */
  do { 
    fgets(line, MAXLINE, fp);
    linenb++;
  } while (line[0] == '%' && !feof(fp));
  WARNING((sscanf(line,"%d %d",&npro,&nedges) != 2,
	   "LIBN_Read code 1 : cannot read the number of processors and number of links, line %d\n",
	   linenb),
	  NULL);
  g->nedges = nedges = nedges*2;

  /* read second line (application ratio) */
  do { 
    fgets(line, MAXLINE, fp);
    linenb++;
  } while (line[0] == '%' && !feof(fp));
  WARNING((sscanf(line,"%d %d %lf",&r_from,&r_to,&r_val) != 3,
	   "LIBN_Read code 2 : cannot read application ratio line %d\n",linenb),
	  NULL);

  /*****************************************/
  /* Init                                  */
  /*****************************************/
  g->npro = npro;
  
  MEM_Needxy(mem, 7*npro, 7*npro,        int   );
  MEM_Needxy(mem, npro, npro,            double);
  MEM_Need  (mem, 4*npro,                double);
  MEM_Need  (mem, 3*npro + 3*nedges + 1, int);

  MEM_Allocate(mem);
  
  realnetcost = g->realnetcost = MEM_Requestxy(mem, npro, npro, int   );
  netcost     = g->netcost     = MEM_Requestxy(mem, npro, npro, int   );  
  g->proratio                  = MEM_Request  (mem, npro,       double);
  proportion  = g->proportion  = MEM_Request  (mem, npro,       double);
  prospeed    = g->prospeed    = MEM_Request  (mem, npro,       double);
  ratio       = g->ratio       = MEM_Requestxy(mem, npro, npro, double);
  
  adjncy      = g->adjncy      = MEM_Request  (mem, nedges,     int   );
  adjwgt      = g->adjwgt      = MEM_Request  (mem, nedges,     int   );
  g->rpro                      = MEM_Request  (mem, npro,       int   );
  g->exec                      = MEM_Request  (mem, npro,       double);
  adjidx      = g->adjidx      = MEM_Request  (mem, npro+1,     int   );
  
  connect     = g->connect     = MEM_Requestxy(mem, npro, npro, int   );
  adjlcy      = g->adjlcy      = MEM_Request  (mem, nedges,     int   );
  latency     = g->latency     = MEM_Requestxy(mem, npro, npro, int   );
  latency1    = g->latency1    = MEM_Requestxy(mem, npro, npro, int   );
  convtxs     = g->convtxs     = MEM_Requestxy(mem, npro, npro, int   );
  convtxs2    = g->convtxs2    = MEM_Requestxy(mem, npro, npro, int   );
  conproc     = g->conproc     = MEM_Request  (mem, npro,       int   );
  
  g->app      = LIBG_Init();

  *adjidx = 0;
  for(i=0;i<npro;i++)
    g->exec[i] = 0;

  /*************/
  /* read file */
  /*************/
  for(i=0;i<npro;i++){
    do { 
      fgets(line, MAXLINE, fp);
      linenb++;
    } while (line[0] == '%' && !feof(fp));
    WARNING((strlen(line) >= MAXLINE,
	     "LIBN_Read code 3 : line too long, line %d\n",
	     linenb),
	    NULL);

    /* proportion */
    WARNING(((tok1 = strtok(line," ")) == NULL,
	     "LIBN_Read code 4 : cannot read processor proportion, line %d\n",
	     linenb),
	    NULL);
    proportion[i] = atoi(tok1);
    if(proportion[i] > maxproportion)
      maxproportion = proportion[i];
    
    /* connectivity  matrix */
    j = *(adjidx++);
    *adjidx = j;

    while((tok1 = strtok(NULL," ")) != NULL && (tok2 = strtok(NULL," ")) != NULL && (tok3 = strtok(NULL, " ")) != NULL){
      j = atoi(tok1);
      dist = atoi(tok2);
      ltcy = atoi(tok3);

      (*adjidx)++;
      *(adjncy++) = j;
      *(adjwgt++) = dist;
      *(adjlcy++) = (int)(ltcy/NBytes);
    }
  }

  /* scale proportion to have fastest processor proportion = 1 */
  for(i=0;i<npro;i++){
    proportion[i] /= maxproportion;
    prospeed[i]   =  1/proportion[i];
  }

  adjidx    = g->adjidx;
  adjncy    = g->adjncy;
  adjwgt    = g->adjwgt;

  adjlcy = g->adjlcy;

  /************************/
  /* generate realnetcost */
  /************************/
  for(i = 0; i < npro; i++)
    for(j = 0; j < npro; j++){
      connect[i][j] = (i==j ? 0 : INFINITY);
      convtxs[i][j] = -1; 
      latency[i][j] = (i==j ? 0 : INFINITY);
      realnetcost[i][j] = (i==j ? 0 : INFINITY);
    }
  for(k = 0; k < npro; k++){
    istart = adjidx[k];
    iend = adjidx[k+1];
    for(i = istart; i < iend; i++){
      x = adjncy[i];
      latency[k][x] = adjlcy[i];
      latency[x][k] = adjlcy[i];
      realnetcost[x][k] = adjwgt[i];
      realnetcost[k][x] = adjwgt[i];
      connect[k][x] = 1;
      connect[x][k] = 1;
    }
  }
  for(k = 0; k < npro; k++)
    for(i = 0; i < npro; i++)
      for(j = 0; j < npro; j++)
        if((connect[i][k] + connect[k][j]) < connect[i][j]){
          connect[i][j] = connect[i][k] + connect[k][j];
          convtxs[i][j] = k;
        }
  for(i = 0; i < npro; i++)
    for(j = 0; j < npro; j++){
      if(convtxs[i][j] == -1)
        continue;
      k = convtxs[i][j];
      if(realnetcost[i][j] == INFINITY){
        a = 0;
        b = 0;
        a = getBandwidth(convtxs, realnetcost, i, k);
        b = getBandwidth(convtxs, realnetcost, k, j);
        if(a >= b){
          realnetcost[i][j] = a;
          realnetcost[j][i] = a;
        }else{
          realnetcost[i][j] = b;
          realnetcost[j][i] = b;
        }
      }
      if(latency[i][j] == INFINITY){
        a = 0;
        b = 0;
        a = getLatency(convtxs, latency, i, k);
        b = getLatency(convtxs, latency, k, j);
        latency[i][j] = a + b;
        latency[j][i] = a + b;
      }
    }
  /**********************/
  /*  generate netcost  */
  /**********************/
  for(i=0;i<npro;i++){
    for(j=0;j<npro;j++){
      netcost[i][j] = realnetcost[i][j];
    }
  }  
  
  /*************************/
  /* generate ratio matrix */
  /*************************/
  g->r_from = r_from;
  g->r_to   = r_to;
  g->r_val  = r_val;
  r_ref = r_val * prospeed[r_from] / netcost[r_from][r_to];
  for(i=0;i<npro;i++)for(j=0;j<npro;j++){
    ratio[i][j] = r_ref * netcost[i][j];
  }
 
  for(i = 0; i < npro; i++){
    for(j = 0; j < npro; j++){
      latency1[i][j] = (int)(r_ref * latency[i][j]);
      if(latency1[i][j] > max_latency){
        max_latency = latency1[i][j];
      }
    }
  }

  if(max_latency > LW_CEIL){
    range = (int)(max_latency / LW_CEIL);
  }

  for(i = 0; i < npro; i++){
    for(j = 0; j < npro; j++){
      latency[i][j] = (int)(latency1[i][j] /range);
    }
  }
  return g;
}

/** 
 * \compute latency weight between processors from and to
 *
 * \param con path recorder 
 * \param netcost bandwidth recorder
 * \param from processor from
 * \param to processor to
 */
int getBandwidth(int** con, int** netcost, int from, int to){
  int k;
  int a, b;
  if(con[from][to] == -1){
    return netcost[from][to];
  }else{
    a = 0;
    b = 0;
    k = con[from][to];
    a = getBandwidth(con, netcost, from, k);
    b = getBandwidth(con, netcost, k, to);
    if(a >= b)
      return a;
    else
      return b;
  }
}

/**
 * \compute latency weight between processors from and to
 *
 * \param con path recorder 
 * \param latency latency recorder
 * \param from processor from
 * \param to processor to
 */
int getLatency(int** con, int** latency, int from, int to){
  int k;
  int a, b;
  if(con[from][to] == -1){
    return latency[from][to];
  }else{
    a = 0;
    b = 0;
    k = con[from][to];
    a = getLatency(con, latency, from, k);
    b = getLatency(con, latency, k, to);
    return a+b;
  }
}

/**
 * \brief compute processor's penalty weights
 *
 * \param graph application graph
 * \param grid computational grid
 */
void LIBN_Progen(GraphType *graph, GridType *grid){
  int i;
  double f, sum0, tmp;
  int totalwgt;

  grid->penalty = U_malloc(grid->npro,int,"LIBN_Progen error code 1\n");
  
  f=0;
  for (i=0; i<grid->npro; i++){
    if (grid->proportion[i]!=0)
      f+= grid->proportion[i];
  }
  f= f/((double)grid->npro);
  sum0 = ((double)graph->nvtxs)/(((double)grid->npro)*f);

  totalwgt = graph->nvtxs;
  
  for (i=0; i<grid->npro; i++){
    tmp = (1-grid->proportion[i])*sum0;
    if (tmp !=0) 
      grid->penalty[i]=(int)floor(tmp);
    else
      grid->penalty[i] = 0;
    totalwgt += grid->penalty[i];
  }
}

/**
 * \brief free the grid structure
 *
 * \param grid structure to be freed
 */
void LIBN_Free(GridType *grid){
  if(grid!=NULL){
    if(grid->penalty)
      U_free("LIBN_Free code 1\n",&(grid->penalty),LEND);
    LIBG_Free(grid->app);
    MEM_Free(&(grid->mem));
    U_free("LIBN_Free code 2\n",&grid,LEND);
  }
}

/**
 * \Initialize the count of number of edges between each pair of processors
 *
 * \param grid computational grid
 * \param graph application graph
 */
void Init_Connect(GraphType *graph, GridType *grid){
  int** connect = grid->connect;
  int **convtxs = grid->convtxs;
  int **convtxs2 = grid->convtxs2;
  int  *conproc = grid->conproc;
  int      npro = grid->npro;

  int*   adjncy = graph->adjncy;
  int*   adjidx = graph->adjidx;
  int*      pro = graph->pro;
  int     nvtxs = graph->nvtxs;

  int i, istart, iend, prov, v, prow, w, j;

  for(i = 0; i < npro; i++)
    for(j = 0; j < npro; j++){
      connect[i][j] = 0;
      convtxs[i][j] = 0;
      convtxs2[i][j] = 0;
    }

  for(v = 0; v < nvtxs; v++){
    prov = pro[v];
    for(i = 0; i < npro; i++)
      conproc[i] = 0;
    istart = adjidx[v];
    iend = adjidx[v+1];
    for(i = istart; i < iend; i++){
      w = adjncy[i];
      prow = pro[w];
      if(prov != prow){
        connect[prov][prow]++;
        if(conproc[prow] == 0){
          convtxs[prov][prow]++;
          convtxs2[prov][prow]++;
          conproc[prow] = 1;
        }
      }
    }
  }
}
